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ABSTRACT 


A  stable  numerical  technique  is  given  for  the  integration  of 
the  fluid  flow  equations  involving  finite  rate  chemical  klntfV.cs. 
The  procedure  permits  a  significantly  larger  integration  step 
than  standard  procedures,  such  as  Runge-Kutta.  Consequently,  the 
amount  of  digital  computer  time  to  numerically  integrate  the 
equations  over  a  specified  domain  is  greatly  reduce!.  To  simplify 
the  implementation  of  the  numerical  procedure  for  generic  chem¬ 
ical  systems,  a  programming  system  has  been  developed  that  auto¬ 
matically  generates  the  necessary  subroutines  from  a  few  inputs. 
The  extension  and  application  of  the  numerical  method  to  partial 
differential  equations  is  also  discussed. 
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1.  INTRODUCTION 

For  the  analysis  and  study  of  flow  problems  involving  finite 
rate  chemical  reactions,  the  equations  of  mass  conservation, 
momentum,  state,  energy,  and  species  conservation  must  be  solved. 
The  latter  equation  although  arpearing  simple  in  form  introduces 
severe  difficulties  in  the  numerical  solution  of  such  problems. 

In  the  case  of  one-dimensional  flow  problems,  the  species  conser¬ 
vation  relationships  for  r.  spcrios  are  represented  by  n-coupled 
non-linear  ordinary  differential  equations  (ODE) .  Numerous 
standard  numerical  methods  can  be  used  to  integrate  these  equa¬ 
tions;  these  procedures  will  be  severely  limited  by  the  size  of 
the  integration  step.  For  two-dimensional  problems,  a  coupled 
system  of  partial  differential  equations  represents  the  species 
conservation  relationship,  and  the  numerical  integration  of  these 
equations  is  likewise  very  difficult  using  standard  explicit  or 
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implicit  methods.  Usually,  the  implicit  method  is  considered 
unconditionally  stable,  and  therefore  an  integration  step  of  any 
size  can  be  used  within  the  limitation  imposed  by  truncation 
error.  Unfortunately,  this  property  of  the  implicit  method 
is  dependent  upon  the  form  of  the  differential  equations. 

Because  of  the  coupling  of  the  n  species  in  each  conservation 
equation,  the  amplification  matrix  is  not  necessarily  convergent, 
and  the  method  can  become  unstable  unless  very  small  steps  are 
used  (see  Section  5). 

Methods  of  overcoming  the  limitation  on  stepsize  have  been 
attempted  by  numerous  investigators.  In  Ref.  7,  the  authors 
develop  a  very  useful  method  of  integrating  a  first  order  differ¬ 
ential  equation  of  chemical  kinetics.  The  procedure  is  limited 
to  uncoupled  equations, and  quite  naturally  other  authors  have 
attempted  to  extend  this  procedure  to  coupled  systems  through 
the  use  of  suitable  transformation!!.  However,  the  transformation 
required  considerable  computation  which  had  to  be  repeated  at 
each  step.  Other  investigators  (Ref.  8)  modified  the  Runge-Kutta 
procedure  and  improved  its  stability  behavior  under  certain 
circumstances.  In  Ref.  9,  a  review  of  some  of  these  approaches 


is  presented. 


•  •  / 


[ 

[ 


For  this  study,  the  method  of  rational  approximations,  which 
was  introduced  by  GASL  (Ref.  3) ,  ha.  been  selected  a.  most  suit- 
sble  for  integrating  the  species  conservation  equations.  The 
details  of  the  method  are  outlined  in  Section  2,  and  it  is  demon¬ 
strated  that  the  method  is  unconditionally  stable.  The  develop¬ 
ment  is  related  to  ODE.  but  in  Section  5  a  method  of  applying 
the  technique  to  partial  differential  equation,  is  discussed. 

The  method  ha.  been  successfully  u.sd  in  numerous  type,  of  prob¬ 
lems  involving  a.  many  a.  31  specie,  and  .'8  reactions. 


es  bel-g  stable,  the  numerical  method  must  be  convenient 
to  use,  if  it  is  to  he  accepted  by  the  engineer  and  scientist. 
Since  the  sethod  of  rational  approximation,  involve,  rather 
lengthy  algebra,  an  algorithm  for  automatically  implementing  the 
«  ti.od  for  a.  generic  chemical  system  ha.  been  developed  a.  part 
of  this  investigation.  Further,  the  algorithm  ha,  been  programmed 
ar.  ,ade  operational  in  conjunction  with  the  GASL  Generalized 
Chemistry  Program,  with  this  sytem  approach,  the  user  merely 
specifies  the  desired  specie,  and  chemical  reactions,  and  the 
necessary  subprograms  are  automatically  generated  by  the  machine. 

A  completely  debugged  subroutine  in  FORTRAN  IV  i„  supplied  as 
output  from  this  system.  The  system  is  very  flexible  and  if 
either  a  specie,  or  a  reaction  is  not  resident  in  the  chemistry 
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library,  it  can  be  readily  added.  The  capacity  of  the  system  is  now 
APO  species  and  200  reactions.  In  Section  3  of  this  report,  the 
details  of  this  system  are  discussed. 

The  advantages  of  using  the  method  of  rational  approximations 
is  demonstrated  in  Section  4;  the  results  for  two  different  disso¬ 
ciating  air  systems  (7  and  11  species)  are  presented.  These  results 
are  compared  with  the  results  from  the  Runge-Kutta  integration 
procedure,  and  the  agreement  is  good.  For  the  problems  studied, 
the  method  of  rational  approximation  required  approximately 
1/20  to  1/25  of  the  machine  time  used  by  the  Runge-Kutta  pro¬ 
cedure.  With  this  very  significant  reduction  in  machine  time, 
the  engineer  can  economically  investigate  large  chemical  systems 
without  resorting  to  approximations  or  simplifications  regarding 
the  number  of  reactions  or  species. 


•  HI 

2.  STABLE  INTEGRATION  PROCEDURES 

A  system  of  n  first  order  ordinary  differential  equations 
can  be  expressed  in  general  form  as 

y'(t)  *  f(t,y)  ;  y(tQ)  -  yQ  (2.1) 

where  y  is  a  vector  valued  function  of  t,  and  f  is  a  vector 

valued  function  of  t  and  y.  Given  the  initial  vector  y  ,  the 

o 

system  of  ordinary  differential  equations  is  to  be  numerically 
integrated  over  the  time  domain  (t^t^). 

In  this  section,  the  numerical  procedures  that  can  be  used 
to  integrate  (2.1)  are  considered.  In  particular  the  stable 
integration  procedures  are  discussed,  j.  or  the  purpose  of  this 
report,  a  numerical  method  of  solution  is  termed  stable  if  an 
error,  once  introduced,  decreases  from  step  to  step  when  the 
eigehvalues  of  the  matrix  [of/dy]  have  negative  real  parts 
(see  Ref.  3) . 

2.1  Discussion  of  Error  Propagation 

Many  methods  have  been  devised  for  obtaining  a  numerical 
solution  of  (2.1).  However,  most  of  these  methods  will  not  give 
a  valid  solution  if  the  stepsize  h  becomes  too  large.  Although 
the  truncation  error  increases  in  a  smooth  manner  as  h  increases. 


the  numerical  solution  becomes  unstable  when  h  exceels  a  certain 
critical  value.  Even  if  the  truncation  or  round-off  errors  are 
small,  the  propagated  error  becomes  unbounded  because  the  ampli¬ 
fication  term  in  the  error  expression  has  a  spectral  radius 
greater  than  one. 

The  unstable  behavior  can  be  illustrated  and  under¬ 
stood  by  considering  the  error  behavior  of  the  Euler  method. 

Lot 


**+?.  -  Yk  +  h\ 


(2.2) 


whero  fk  -  fft^,  yk)  h  is  the  stepsize  between  tfc+1  and  t  ,  and 
Sc  ■  kh.  The  stepsize,  h,  can  be  assumed  a  constant  for  this 
discussion.  The  true  solution  can  be  expressed  as 


k+1 


yk  +  hfk  +  Vi 


(2.3) 


where  rk+1  is  the  truncation  error.  By  subtracting  (2.2)  from 
(2.3),  applying  the  mean  value  theorem,  and  letting  e 
the  error  expression  becomes 


k  "  *k  ’  V 


k+1 


(I  +  hr 

l  ay 


+  r. 


k+1 


(2.4) 


df 


where  ~  is  the  matrix 


iX  l&yj]  *  The  operator  on  the  error  terr. 
%  is  the  amplification  matrix  which  determines  the  stability 

behavior  of  the  system  of  liquations.  Dropping  the  truncation 
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<A  * 


error  and  assuming  that  is  a  constant!  an  initial  error  vector 


by 


eQ  will  propagate  as 


'k+1 


-(-•»«) 


k+1 


(2.5) 


If  the  eigenvalues  of  ~~  have  positive  real  parts,  the 
solution  would  be  increasing  and  the  numerical  error  cculd  like¬ 
wise  be  expected  to  increase.  If,  however,  the  eigenvalues  all 
have  negative  real  parts  and  an  error  increases  from  step  to 
step,  the  numerical  solution  will  then  be  partially  unstable. 

If  the  error  generated  by  Euler’ a  method  is  not  to 

increase,  the  eigenvalues  of  the  amplification  matrix  I  +  h  — 

by 

must  be  less  than  one  in  absolute  value.  Otherwise,  as  the  matrix 
is  raised  to  higher  powers,  the  norm  of  the  error  vector  will 
continually  grow.  Let  the  complex  eigenvalue  of  be 


"  +  Lfil 


where  a^  >  0.  Then  by  imposing  the  stability  requirement,  a  suf¬ 
ficient  condition  for  stability  is 


h  <  2/a 


i 


8 


0k/ 


If!  ot  is  large,  as  in  the  chemical  kinetic  case,  this 
requirement  may  result  in  a  small  h  and  a  large  amount  of  computing 
time  if  the  domain  of  integration  is  large.  The  above  behavior  for 
the  Euler  method  is  characteristic  of  most  standard  single  and 
multi-step  methods.  Extending  the  analysis  to  the  more  complica¬ 
ted  methods  does  not  present  any  difficulty  other  than  cumbersome 


algebra  and  has  been  reported  in  Ref.  4.  For  the  purpose  of  later 


discussion,  the  stability  criterion  for  the  fourth  order  Rung- 
Kutta  method  is 


h  <  2,8 


where  X  is  an  eigenvalue  of  maximum  modulus. 

The  equations  describing  the  chemical  reactions  for 
dissociating  air  can  be  used  to  illustrate  the  restrictions  on 
stepsize.  For  a  streamtube  type  of  calculation  for  a  reentry  body 
(see  Section  4),  the  dominant  eigenvalue  is  about  vxlO4  ,  and  for 
the  Euler  method  to  remain  stable,  h  <  ,6  x  10-4  (ft).  To  inte¬ 
grate  the  equations  ov»r  a  distance  of  0.5  ft  would  require  about 
8000  steps.  These  small  steps  are  required  even  though  the  solu¬ 
tions  change  very  little  from  one  step  to  the  next. 


machine  time)  would  be  required.  As  shown  in  the  next  subsection 


the  limitation  can  be  removed  by  using  numerical  methods  which 


are  unconditionally  stable 


2.2  The  Pade'  integration  Method 


A  class  of  stable  integration  formulas  can  be  derived 


by  using  known  properties  of  linear  ODE.  These  formulas,  while 
directly  applicable  to  a  system  linear  ODE,  may  be  applied 
to  systems  of  non-linear  differential  equations  by  linearizing 
at  each  step.  Even  considering  the  additional  computer  time 


necessary  to  linearize  the  equations,  the  total  time  needed  to 


solve  a  complete  problem  is  mjch  less  than  required  by  any 


standard  method 


Consider  a  system  of  n  linear  ODE  with  constant  coef 


The  formal  solution  of  this  oquatior  at  t*h  is 


A. 


» 


*  t 


In  Eq.  (2.7),  the  exponential  term  can  be  appro:. imated 
in  numerous  ways,  and  the  form  of  the  approximation  will  determine 
whether  the  procedure  is  stable.  For  our  purposes,  a  convenient 
method  of  approximating  any  function  is  the  rational  approxima¬ 
tion,  which  can  be  derived  using  the  Pade'  transformation  of  the 
power  series.  If  P  and  Q  are  polynomials  of  degree  p  and  q, 
respectively,  the  exponential  is  represented  as 


’iA  —i 

e  -  Q  P  +  E  (hA) 


The  rational  approximation  agrees  with  the  power  series  of 
exp  (hA)  for  at  least  p  +  q  +  1  terms  and  the  residual  error  is 
denoted  by  E (hA) .  The  polynomials  P  and  Q  for  the  exponential 
function  (Ref.  5)  are  given  by 


r  (p  +  g  -  k) !  p! _ 

k-0  *P  +  '  k*  ' 


(hA) 


(P  +  q  -  k)  J  q.' 


(p  +  q)  I  k!  (q  -  k) ! 


7  (-hA) 


Substituting  the  rational  approximation  into  formula 
(2.7)  yields  the  desired  integration  formula: 


y (h)  “  Q-  P  (yQ  +  A  b)  -  A-1  b 


Or  rearranging  the  terras. 


Equation  (2.8)  ircludee  numerous  types  of  single  step 
integration  methods  and  not  all  of  these  methods  are  stable. 

The  Euler  method,  which  has  been  shown  to  be  partially  stable 
in  the  previous  section,  is  a  special  case  of  (2.8)  with  Q  =  I 
and  P  ■  I  +  hA  (q  ■  0  and  p  ■  1) .  In  this  report,  only  the 
stable  procedures  are  of  interest  and  for  these  methods  q  *  p. 

Of  particular  interest,  are  those  methods  where  the  diagonal 
(P  ■  q)  of  the  Pade'  table  is  used.  Such  methods  shall  be  called 
Pade'  integration  procedures  or.  rational  approximation  methods 
in  this  report.  For  the  cases  p  *  q  -  1  and  p  ■  q  ■  2,  the 
following  expressions  are  obtained  from  (2.e): 

y(h)  -  (I  -  ~  hA)"X  [  (I  +  “  hA)yQ  +  hb]  (2.9) 

y(h)  -  (i  -  j  hA  +  y-  (hA) 8 1  1  |jl  +  y  hA  +  yy  (hA)a)yQ  +  hbj 

Since  these  formulas  must  agree  with  the  power  series 
for  at  least  p  +  q  +  1  terms,  their  truncation  errors  are 
immediately  seen  to  be  0(ha)  and  0(hB),  respectively.  Similarly, 
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integration  formulas  can  be  derived  for  any  order  of  truncation 
error. 

The  stability  of  the  Pade'  integration  method  can  be 
investigated  by  considering  the  following  form  of  (2.8) 


y  (h)  -  Q-1  PyQ  +  Q-1  Rb 


where  R  represents  some  polynomial  in  (hA) .  Applying  the  method 
in  a  step  by  step  manner  yields 


Q_i  P  yk  +  Q"1  Rb 


ffc+1  ~  w  *  Xy  T  w  (2.10) 

where  y^  *  y(kh) .  The  true  solution  y  can  also  be  expressed 


as 


rk+l 


o’1  Rb  ♦  Tk+1 


(2.11) 


where  T^+1  is  the  truncation  error.  Denote  the  error  at  each 
step  by 


®k  "  yk  "  *k 

and  subtract  (2.10)  from  (2.11)  to  obtain 


-i 


Vi m  (a  p>\  +  T 


k+1 


(2.12) 


In  order  to  see  how  an  error  from  any  source  propagates, 
we  may  neglect  the  truncation  error  term  and  assume  that  a  single 


0 

0 

D 

0 

D 

0 

D 

0 

D 
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error  eQ  is  conmittad  at  one  step,  say  the  first  step.  Equation 
(2.12)  then  becomes 

®k  "  Q  1  P  ek-l  "  (Q  1  P)3  ®k-2  ’  "  (Q-1  p)k  ®Q  (2.13) 

After  k  steps  an  initial  error  eQ  has  increased  or  decreased, depending 

upon  the  amplification  matrix  (Q  P) . 

From  the  definition  of  stability,  when  the  eigenvalues 

of  the  matrix  A  have  negative  real  parts,  e^  must  decrease  as  k 

increases.  Since  e^  is  a  constant,  e^  can  decrease  only  if  the 
— i  )c 

matrix  (Q  P)  is  convergent;  that  is,  the  eigenvalues  of 
(Q  1  P)  must  be  less  than  one  in  absolute  value.  For  every 
eigenvalue  X^  of  the  matrix  A  in. (2.6),  there  is  a  corresponding 
eigenvalue  for  Q_1  P.  For  convenience,  the  stability  criteria 
‘is  stated  in  terras  of  X  by 

i  , 

|!p(X)/q(X)|  <  l 

.  .  i  '  1 

Thus  the  stability  requirement  becomes 

.  .  .  •  ,  i 

|P(X)|  <  |Q(X)|  (2.14) 

where  the  complex  eigenvalue  X  -  hX  ■  -  o  +  ifi  (a  >  0) .  For  p  *  q 
it  can  be  shewn  that  the  stability  requirement  is  satisfied. 


+  *0 


l 


%  ' 


The  Pade '  integration  method  ia  directly  applicable  to 


linear  systems  of  ordinary  differential  equations.  Some  additional 


steps  are  required  to  apply  the  method  to  chemical  kinetic  equa¬ 


tions,  which  are  represented  by  (2.1).  The  non-linear  system 


of  equations  can  be  reduced  to  the  form  (2.6)  by  expanding  the 


function  f  in  a  Taylor's  Series  and  neglecting  high  order  terms, 


y  "  f(to'yo)  (to'yo>  *y  +  0(Ay3> 


(2.15) 


where  by  -  y  -  yQ  and  ^  (tQ,yo)  denotea  the  matrix  A  in  (2.6) 


with  elements 


ij  ay. 


Thie  indices  i  and  j  denote  the  ifc^  and  elements  in  the  f  and 


y  vectors,  respectively.  Hence  grouping  terms  as  in  (2.6) 


y  '  9y  (to'yo,y  +  (f(to'yo>  -  a?fro'yo,yo)  +  0(Ayi,> 


Ay  +  b  +  0(Aya) 


(2.16) 


Since  the  equation  must  be  linearized  at  each  step  of  the  inte¬ 


gration,  the  residual  error,  which  is  denoted  by  ^+1  *  0(Ay3)  for 


the  k  step,  must  be  introduced  in  Eq.  (2.12), 
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Vi  -  (Q'1  p,ek  +  0-1  R  \+i  +  Tk+i 


(2.17) 


Again,  assuming  the  matrix  A  and  the  vector  b  are  nearly  constant 
over  (k+l)  steps  then 


.  -i  k+l 
*k+i  *  (Q  p) 


k+l 


,so  +  l  (Q'‘  P)’1'  <Q"‘  Rlu  +  V 

v-1 


(2.18) 


Equation  (2.18)  which  is  similar  to  (2.13)  includes  the  influence 

of  the  truncation  error  and  the  residual  error  of  linearization. 

Again,  the  amplification  matrix  Q  1  P  must  be  a  convergent  matrix 

**i  k 

for  a  stable  system.  The  rate  of  decay  of  (Q  P)  and  the 
magnitude  of  the  error  Q  1  R  introduced  at  each  step  will 

determine  whether  efc+^  is  decreasing  with  k.  Since  the  linear¬ 
ization  and  the  truncation  errors  appear  together,  the  integration 
method  for  a  particular  type  of  problem  should  be  selected  to 
achieve  a  truncation  error  approximately  the  same  size  as  the 
linearization  error.  Hence,  if  a  large  residual  error  is  present, 
a  low  order  integration  formula  should  be  used. 


/. 
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3. 


IMPLEMENTATION  OF  THE  INTEGRATION  PROCEDURES 
As  shown  in  Section  2,  the  Pade'  approximation  technique  for 
solving  systems  of  ordinary  differential  equations  is  directly 
applicable  to  linear  equations.  For  non-linear  systems  of  equa¬ 
tions,  the  functions  f.(yj(t),i(  j  «  1,  ...  n  where  n  is  the 
number  of  species  in  the  chemical  system  must  be  expanded  in  a 


Taylor's  Series.  The  terms  of  second  order  and  higher  are 
neglected  in  the  series  for  each  function.  At  each  step  of  the 
integration,  this  linearized  form  of  the  functions  must  be 
evaluated.  Hence,  when  the  Pade'  method  of  integration  is 
implemented  on  a  digital  computer  it  is  convenient  to  develop 
two  separate  routines.  The  first  routine  calculates  the  linear¬ 
ized  form  of  the  functions  and  the  second  routine  integrates 
the  linear  form  over  a  single  step  in  the  independent  variable. 


Using  the  terminology  of  Section  2  [Eq.  (2.16)1,  the  former 
routine  calculates  the  matrix,  A,  and  forcing  vector,  b,  whereas 
the  latter  routine  solves  a  system  of  linear  first  order  ordinarv 
differential  equations. 

The  system  of  linear  ordinary  differential  equations  is 

defined  by  the  matrix,  A.  The  routine  which  solves  this  system 

, 

is  related  to  the  chemistry  problem  only  through  the  size  of  the 
A  matrix.  Hence,  this  integration  routine  need  only  be  written 


purpose  integration  routine  has  been  written  and  combined  with 
the  one -dimensional  flow  equations  to  form  a  general  program  for 
the  analysis  of  strearatube  problems  involving  chemical  kinetics. 
This  program  and  its  use  is  described  in  Subsection  3.3. 

The  subroutine  for  forming  the  matrix,  A,  and  forcing 
vector,  b,  is  directly  dependent  upon  the  types  of  reactions  and 
species,  and  consequently,  this  routine  must  be  recoded  for  each 
new  chemistry  system.  For  a  large  number  of  reactions  and 
species,  the  analysis  and  coding  ace  very  tedious  and  error  prone. 


To  overcome  this  difficulty  an  algorithm  has  beer,  developed  for 
forming  A  and  b  given  a  generic  chemistry  system.  Furthermore, 
the  algorithm  was  implemented  in  conjunction  with  the  Generalized 
Chemistry  Program  (GCP)  developed  at  GASL  (see  Ref.  1).  Given  a 
small  number  of  inputs,  the  GCP  generates  in  FORTRAN  language 
the  subroutines  which  can  be  used  with  the  integration  program. 

In  Section  3.1,  the  algorithm  for  forming  the  A-raatrix  and 
b-vector  is  described  along  with  the  programming  procedures.  A 
general  descrip  Ion  of  the  GCP  is  given  in  Section  3. 2, and  the 
required  input  car  Is  for  generating  the  subroutine  are  described. 
The  GCP  has  been  successfully  used  on  a  variety  of  cheraibtry 
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problems  involving  as  many  as  31  species  and  78  reactions 


3.1  Algorithm  for  the  Linearized  Form 


Problems  in  chemical  kinetics  give  rise  to  differential 


e 


equations  of  the  form: 
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i«l, 2, . . .n  (3.1) 


where  n  is  the  number  r£  reacting  species,  m  is  the  total  number 


of  reactions  (the  sum  of  forward  and  reverse  reactions),  is 


the  number  of  moles  per  unit  mass  of  the  mixture,  and  the  dot 


denotes  differentiation  with  respect  to  time.  T  and  p  are  the 


temperature  and  density,  respectively.  The  stoichiometric 


coefficients  and  are  defined  by  the  following  expression 


for  the  j  chemical  reaction  equation: 


i  Ki  \  - 1  pi 


(3.2) 


where  and  are  the  reactants  and  products,  respectively. 


The  reao-ion  rate  for  the  j  reaction  is  k^  which  is  taken  as 


a  function  of  temperature  only  and  has  the  form 


k_,  (T)  *  a .  T  exp  (-e./T) 
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where  a  ^ ,  b_, ,  and  e^  are  constant*. 


The  right  side  of  Eq.  (3.1)  represents  the  function  f 


which  must  be  expanded  in  a  Taylor's  Series  as  shown  in  Eq.  (2.15) 


In  terms  of  the  functions  fi  the  elements  of  the  matrix,  A,  and 


vector,  b,  are 


is  3y 


(3.3) 


bi  ■  fi(to'yo>  -  l  aiit<yo): 

k-1 


(3.4) 


where  index  i  denotes  the  species  production  equation,  index 


denotes  the  species,  and  (yQ)^  is  an  element  of  the  vector  of 
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species  at  time  tQ.  The  algorithm  which  gene 
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for  (3.3)  and  (3.4)  requires  numerous  constants  (e.g.  a^,  b^ 


e j  •  v±y  etr-’)  •  of  the  required  constants  are  on  the 


library  tape  of  the  GCP  which  was  developed  by  GASL.  The 


algorithm  for  handling  the  constants,  defining  the  arrays,  and 
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operating  upon  these  arrays  is  described  below. 


When  a  particular  set  of  > pecies  is  selected,  tne 


species  are  consecutively  numbered  accrrding  to  their  order  of 


input.  Each  speciesis  then  repissitnl  d  ’'v  the  corresponding 


integer  throughout  the  logic  of  the  generator  program.  These 
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integers  are  used  in  a  12  x  m/2  array,  which  represents  the  reac¬ 
tions  (3.2).  The  first  six  columns  of  the  a,“ray  are  used  for 
the  left-hand  side  of  the  1  sections  and  the  last  six  columns 
represent  the  right-hand  side.  The  integers  corresponding  to 
species  with  non-zero  stoicniometric  coefficients  in  the  reaction 
equation  are  placed  in  these  columns.  If  a  stoichiometric  coef- 
f icient  is  greater  than  one,  the  specie  number  is  simply  repeated 
the  proper  number  of  times.  For  example,  the  reaction 

F«  +  2Ri,  +  Ri.  -  P7  +  Pa8 

would  occupy  one  row  of  the  reaction  array  and  would  nppear  in 
storage  as  (each  box  is  a  storage  cell) s 


5 

12 

12 

18 

Note  the  reactant  a  occupies  two  memory  cells  because  its 
stoichiometric  coefficient  is  two  in  the  reaction  equation.  Sin 
all  reactions  are  considered  reversible  (if  not  the  reverse  reac 
tion  rate  is  merely  set  equal  to  zero) ,  the  reverse  reaction  is 
obtained  by  merely  scanning  the  row  from  right  to  left. 

As  indicated  above,  the  current  form  of  the  logic  of 
the  generator  program  is  restricted  to  reaction  equations 
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involving  6  or  leaa  reactants.  Four  of  these  reactants  may  be 


different  species.  The  aame  restrictiona  are  placed  upon  the 


products.  Withia these  restrictions  our  survey  indicates  that  all 


reactions  of  current  interest  in  the  area  of  reentry  physics 


can  be  included.  If  in  the  future  more  sophisticated  chemical 


reactions  are  required  (e.g.  5  body  reactions  or  reactions 


involving  large  stoichiometric  coefficients) ,  the  program  logic 


is  readily  extended. 


The  reaction  array  ia  subdivided  into  three  distinct 


groups.  Each  reaction  which  is  to  be  stored  in  memory  is  first 


classified,  and  then  assigned  to  the  corresponding  group.  The 


first  group  consists  of  all  of  the  noncatalytic  reactions.  The 
second  group  is  for  nondesignated  catalytic  reactions;  that  is. 


the  reaction  rate  ia  independent  of  the  species  which  are  used 


as  catalysts.  lu  the  third  group,  the  designated  ca^-alytic 


reactions  appear.  For  these  reactions,  the  reaction  rate  ia 


dependent  upon  the  species  designated  as  a  catalyst.  There  may 


be  p(p  £  10)  different  reaction  rates  for  each  designated 


catalytic  reaction.  Consequently,  the  set  of  species  must  be 


partitioned  into  p  or  leas  subsets  for  each  reaction  in  the  third 


group.  An  auxiliary  array  is  used  co  define  these  subsets  and 


their  corresponding  reaction  rates. 
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All  the  information  that  is  necessary  to  form  the  func* 

tions  f^  and  then  differentiate  them  with  respect  to  each  of  the 

species  as  implied  by  Eqs.  (3.3)  and  (3.4),  is  contained  in  the 

arrays  described  above.  These  arrays  are  initially  formed  by 

the  program  generator,  and  then  processed  as  shown  in  the  flow 

diagram  on  page  24  to  encode  the  desired  program.  There  are  two 

* 

variations  of  the  object  program  depending  upon  whether  or  not 
the  derivative  of  f^  with  respect  to  temperature,  T,  is  needed 
(see  Subsection  3.3).  If  they  are  needed, an  array  for  these 
values  is  added  to  the  dimension  cards.  Coding  is  then  generated 
to  define  the  reaction  rate  constants,  calculate  powers  of  p 
needed  in  subsequent  calculations,  and  calculate  the  reaction 
rates  for  the  given  temperatures.  Coding  io  also  generated  to 
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initialize  the  A  array  to  zero, and  if  the  temperature  derivative 


is  needed  an  appropriate  indicator  is  generated. 

The  next  phase  in  the  flow  diagram  is  the  generation 
of  the  code  for  the  partial  derivatives  of  each  reaction  with 
respect  to  each  species  that  has  a  nonzero  stoichiometiic  coef¬ 


ficient.  This  coding  is  used  to  define  the  two-dimensional 
array  PL.  Coding  is  then  produced  to  evaluate  the  terms  of  the 


A  matrix  from  combinations  of  the  terms  in  the  PL  arrav.  Th*» 
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terms  (T)  p  (py^)  are  then  encoded,  and  the  coding  pro¬ 

duced  to  form  the  f  At  this  point  in  the  program,  either  the 
subroutine  is  terminated  with  a  RETURN  statement  or,  if  required 
the  terms  df^/dT  are  encoded  prior  to  the  RETURN. 

The  logic  of  the  generated  object  program  is  shown  in 
the  flow  diagram  on  page  27.  Following  the  initialization 
phase,  the  program  evaluates  the  reaction  rates  for  a  given 
temperature.  This,  however,  is  only  done  if  the  temperature  is 
sufficiently  different  from  the  temperature  at  &  previous  entry. 


Otherwise  the  old  values  are. used.  The  partial  derivatives  of 
-i  n 

the  term  k , (T) p  n  (py.)  with  respect  to  the  y,  with  non- 
j  <t«l  v  t 

zero  v'  are  then  evaluated  and  stored  in  the  PL  array.  These 

terms  are  then  combined  to  form  the  matrix  A.  The  terms 
n  i/' 

-1  P  ci 

(T)  p  fl  (py. )  are  then  formed  and  combined  in  an  array 
J  -6=1  ^ 

F,  which  corresponds  to  f^.  If  the  temperature  derivatives  are 
needed,  the  derivatives  of  these  terms  wich  respect  to  T  are 
evaluated  and  stored  in  the  vector  D. 

A  listing  of  the  program  generator  and  of  ^  typical 
generated  program  may  be  found  in  Appendices  A  and  B. 


I 


•V  „  •  ,  *• 


■-  "•  ■  '  >•,  4*  r 

N.  'V  3  ' 


FLOW  DIAGRAM  FOR  THE  GENERATOR  PROGRAM 


INITIALIZE 


Count  Number  of  Seta 
of  Reaction  Rate  Data 


Tempera tur 


Derivative 


Generate 
Dimension  Cards 


Generate 
Dimension  Cards 


Generate  Reaction  Rate  Data  from  Library  Tape 


Generate  Coding  for: 


a)  Calculation  of  Powers  of  p 


b)  Calculation  of  k.. (T) 


c)  Zeroing  out  of  A  Array 


Temperature 
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Generate  and  encode  the  following; 


*?k  [kj !T)  p’'  £  «■»*>  <J]  -  ««>*> 


all  species  with  a  non-zero  stoi¬ 
chiometric  coefficient  in  reaction 

<j) 
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Generate  and  encode  the  A  array: 
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Generate  and  encode  the  followings 


w(j)  -  (T^p“l  n  (py.) 
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FLOW  DIAGRAM  FOR  THE  GENERATED  PROGRAM 
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3.2  Use  of  the  General  Chemistry  Program  (GCP) 


In  this  subsection,  the  procedure  for  generating  a  pro¬ 


gram  for  the  linearized  model  of  the  production  equations  for 


finite-rate  chemical  systems  is  described.  Emphasis  is  placed 


on  the  input  format,  operating  instructions  and  documentation  of 


the  GCP.  A  description  of  the  GCP  together  with  sample  inputs 


and  outputs  is  also  included. 


ifrfab 


3.2.1  General  Description  of  the  GCP 


The  GCP  is  a  modular  series  of  programs  which  interact 


to  generate  other  programs.  As  shown  in  Fig.  1,  after  accenting 


the  input  quantities  the  system  can  take  three  general  types  of 


action.  The  supervisor  program  will  select  a  subprogram  to 


modify  or  update  the  library  tape,  to  generate  a  report  or  to 


generate  one  or  more  FORTRAN  IV  subroutine  source  decks.  The 


former  two  options  are  actually  support  functions  for  the  latter 


option?  that  is,  the  generation  of  source  decks  is  the  primary 


objective  of  the  GCP.  Reporting  upon  the  contents  of  the  source 


decks  and  maintaining  a  library  of  chemical  kinetic  data  are 


necessary  activities  to  achieve  this  objective. 


The  library  tape  contains  data  for  various  species  and 


reactions  involving  these  species.  For  each  species,  the 


library  tape  contains  its  chemical  symbol,  molecular  weight 
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cubic  polynomial  fits  of  enthalpy  and  entropy  for  various  tempera¬ 
ture  runges,  and  coefficients  for  mixture  viscosity,  binary  diffus¬ 
ion,  specific  heat  and  thermal  conductivity  models,  for  each 
chemical  reaction,  the  library  stores  the  complete  react' -in  equa¬ 
tion,  the  forward  and  backward  reaction  rates,  and  an  indicator 
that  denotes  if  the  reaction  is  catalytic.  If  the  reaction  xs 
catalytic,  it  is  further  classified  into  Ohfe  of  two  groups.  In 
one  of  the  groups  alx  species  act  as  attbalysts  with  a  common  reaction 
rat j;  in  the  other  group  the  Bpecies  (specified  by  the  user)  are 
separated  into  ten  or  less  designated  sets  and  each  set  has  a 
different  reaction  rate. 
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The  library  can  be  modified  by  providing  the  proper 
control  cards  to  the  supervisor  program  as  indicated  in  Fig.  1. 
Then  the  program  control  is  given  to  the  library*  update  program 
which  reads  the  new  information  or  modifications  of  data  into 
the  primary  memory  of  the  computer.  JL  new  master  library  tape 
is  then  created  by  merging  the  original  tape  and  the  modifica¬ 
tions.  Upon  completion/  the  program  documents  all  the  modifica¬ 
tions  andr  if  requested/  a  complete  report  of  the  library  is 
provided  using  the  report  generation  programs  Cor trol  is  then 
returned  to  the  supervisor  program. 


The  report  generation  program  serves  to  document  the 
contents  of  the  library  (as  indicated  above)  or  the  contents  of 
thci  generated  program.  The  purpose  of  the  documentat ion  is  to 
provide  hard  copy  for  reference  on  a  particular  program  or  for 
the  inclusion  in  a  technical  report.  As  a  convenience,  the  format 
of  the  report  firs  upon  the  standard  size  paper  ( 8*5  x  11  inches) . 

A  typical  output  is  illustrated  by  the  generated  subroutine  FG 
in  Appendix  B. 

The  third  option  of  the  GCP  is  the  generation  of  a 
source  program  in  FORTRAN  IV  language.  As  indicated  in  Fig.  1, 
the  system  can  generate  three  distinct  types  of  programs: 
equilibrium  chemistry,  finite  rate  cnemistry,  and  the  linear¬ 
ized  form  of  the  finite  rate  equations.  In  this  report,  only 
the  latter  option  is  of  prim***;'  interest,  but  for  completeness, 
a  brief  description  of  the  other  two  options  is  included. 

The  equilibrium  chemistry  yvierator  of  the  GCP  is 
used  to  automatically  develop  subro^f  ,*»es  which  will  compute  the 
equilibrium  state  of  a  chemistry  system  given  the  pressure  rnd 
the  temperature.  The  equilibrium  state  is  calculated  by  mini¬ 
mizing  the  Gibbs  Free  Energy  function.  The  minimization  is 
accomplished  by  numerical  iteration  using  the  Newton-Raphson 
method  in  n  variables.  This  method  is  readily  extended  to  the 


pressure-entropy  problems  in  a  direct  way  as  described  in  Ref.  2. 

The  finite  rate  chemistry  option  is  intended  to  gen- 

. 

erate  programs  which  evaluate  in  an  exact  manner  the  right  side 
of  Eq.  (3.1)  for  a  specified  set  of  species  and  reactions.  By 
specifying  a  small  number  of  inputs  (e.g.  the  species  and  reac¬ 
tion  numbers  as  they  appear  on  the  library  tape),  the  system 
will  produce  the  FORTRAN  IV  source  program,  which  will  include 
subroutines  for  computing  the  mixture  enthalpy 
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where  o i  is  the  i  mass  fraction,  and  is  the  polynomial  fit 

til 

of  enthalpy  for  the  i  species.  The  coefficients  of  the 


polynomial  h^  are  taken  directly  from  the  library  tape  of  the 
GCP.  In  addition,  this  portion  of  the  system  also  generates 
a  subroutine  which  will  compute  the  temperature  gi jen  the 
enthalpy.  Numerous  other  functions  can  be  automatically  gener¬ 


ated  under  this  option,  and  these  definitions  and  details  can  be 
found  in  Ref.  1. 

The  third  type  of  program  generated  is  the  linearized 
form  of  Eq.  (3.1).  A  source  program  in  FORTRAN  IV  is  auto¬ 


matically  developed  to  evaluate  the  elements  of  the  A  matrix  as 
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defined  by  Eq.  (3.3).  Further,  a  vector  with  element*  ia  also 


generated,  and  this  vector  and  the  elements  of  the  A  matrix 


can  then  be  used  by  the  calling  program  to  compute  the  vector  b 


in  Eq.  (3.4).  The  method  used  to  generate  the  program  is  des¬ 


cribed  in  Section  3.1  and  samples  of  a  typical  generated  program 


can  be  found  in  Appendix  B. 


3.2.2  Input  Format  for  the  Generation  of  the 
Linearized  Source  Program 


1 


This  section  describes  the  input  format  for  the  linear- 


ized  option  of  the  GCP.  In  the  following  table  of  card  formats, 


the  cards  numbered  1  through  9  are  required  for  the  linearized 


option.  Card  number  3  is  used  to  specify  whether  the  generated 
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routine  is  to  include  the  logic  for  iteration  on  t  ".perature. 


The  cards  10  through  14  are  used  to  encode  a  subroutine  to  com¬ 


pute  : 


„J  MW 

\  * 


(1)  Mixture  enthalpy  given  temperature  and 


v  *  , 


*  -  • 


species  mass  fractions. 


(2)  Temperature  given  mixture  enthalpy  and 


’  * 


species  mass  fraction. 


*  r  •  * 


This  rubroutine,  which  is  not  specifically  a  part  of  the  linear 


.  » 


ized  subroutine,  is  usually  required  by  the  main  program,  and 


for  completeness  the  me the d  of  obtaining  this  routine  is 
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i 
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ard  No. 


4\4# 


. 


1 


Column  i 


41-45 


46-50 


1-72 


. 


Format 


6A12 


( 


_ Description 


month 


date 


year  ) 


ns  -  number  of  species  to  be 
included  £  ICO 


,  . 


nr  «  number  of  reactions  £  100 


< 


■  * 


'0"  if  temperature  iteration  ii 
desired,  else  "l" 


%  • 


Chemical  symbols  of  the  species  to 
be  included  are  punchod  six  species 
per  card.  Each  species  is  allo¬ 
cated  12  columns.  The  symbols  can 
appear  anywhere  *n  the  12  columns. 
The  order  of  the  species  on  these 
cards  determines  the  order  of  the 
species  in  the  program.  Species 
by.nbol  spelling  must  agree  with 
the  library  report. 


» /*  - 

< 


•  / 

■  •  i 


u  • 


v 

/  ’ 


vV 

T 


- « 


:  t 


mm 


•  ■ 


r.  • 


.  r 


V 

•  I 


It* 


v 


* 


Card  No.  ( olumns  Format 


•  r  ^  r  ^ 

>,  5  ,5 


1-70 


41- $5 


13, 13 13' 


1415 


KM 

rary 


Reaction  numbers  from  the  library 
report  of  the  reactions  to  be 
included  are  punched  fourteen 
reactions  per  card.  The  numbers 
are  right  adjusted  and  must  be  in 
ascendina  order. 


5 


Four  blank  card^ 


Same  as  card  #1 


Same  as  columns  41-45,  card  2 


Columns  3,  4,  5,  6,  9,  14,  15,  18, 
19,  21  contain  "1";  all  others  are 
blank. 


Same  as  cards  4,  4',  4' 
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3.2.3  Operating  Instructions  for  Generating 
the  Linearized  Source  Program 


The  set  of  data  card*  described  in  Subsection  3.2.2  is 
placed  between  the  two  end-of-file  cardj  at  the  back  of  the  GCP 
running  deck.  The  running  deck  is  then  executed  on  an  IBM 
7090/94  under  version  13  of  IBSYS.  The  instructions  to  the 
machine  operator  are 3 

(1)  Mount  the  GCP  library  tape  on  A5y  mount 
scratch  tapes  on  A6,  B5,  B6. 

(2)  Execution  time  is  less  than  5  minutes. 


(3)  When  job  terminates,  print  Bl  and  punch 
two  files  from  A6. 


3.2.4  Output  of  the  GCP 

The  output  from  the  GCP  consists  of  a  printed  report 
of  the  program  and  a  source  program  deck  in  FORTRAN  IV.  In 
Appendix.  B,  the  generated  subroutine  for  the  linearized  form  is 
listed  under  the  subprogram  r*G.  This  program  is  for  dissociating 
air  consisting  of  the  following  eleven  species:  0,  Oa  ,  N,  Na , 
N0+,  e”,  NO,  o\  O* ,  N+,  and  Na.  The  subroutine  FT,  which  is 
also  automatically  encoded,  is  used  to  compute  mixture  enthalpy 
or  conversely  temperature  given  mixture  enthalpy. 


'  '  *  * 
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The  inputs  required  to  generate  subroutines  FG  and  FT 
are  listed  at  the  rear  of  Appendix  R  and  can  be  compared  with  the 
description  under  Subsection  3.2.2. 

3.3  Streamtube  Program 

Having  developed  a  computer  program  which  will  auto¬ 
matically  produce  a  routine  that  evaluates  the  linear  approxi¬ 
mation  to  the  differential  equations  of  any  chemical  system,  it 
was  found  desirable  to  write  a  general  program  that  utilizes 
these  routines.  A  subroutine  was  written  that  uses  the  Pade ' 
method  to  integrate  any  system  of  linear  ordinary  differential 
equations.  This  subroutine  is  combined  with  the  linearization 
routine  and  a  routine  to  obtain  the  temperature  from  the  given 
temperature-enthalpy  relationships;  the  combined  program  is 
used  to  calculate  the  chemical  kinetic  properties  along  the 
streamtube.  The  program  can  be  adapted  to  any  chemical  kinetic 
system  by  merely  changing  the  linearization  and  temperature- 
enthalpy  -outines,  which  are  both  obtained  from  the  GCP. 

The  Pade'  integration  routine  used  in  the  program 
assumes  that  y(0)«0.  This  assumption  simplifies  Eqs.  (2.9)  and 
the  linearized  system  is  then  conveniently  expressed  in  terras 


of  by  as 


(I  -  “  hA  +  KA" ) by  =  hb 


(3.5) 


- 


.  f . 


.  *. 


where  K«0  for  a  second  order  integration  scheme  and  K“h’/12  for 

a  fourth  order  integration  echeme.  By  permitting  K  to  be 

specified  as  an  input,  either  the  second  order  or  fourth  order 
« 

methods  can  be  conveniently  selected.  In  cases  where  the  higher 
order  scheme  offers  no  great  advantages,  the  elimination  of  the 
calculation  of  the  matrix  A8  can  save  an  appreciable  amount  of 
computer  time.  The  system  (3.5)  is  solved  by  the  Crout  reduc¬ 


tion  method  to  obtain  Ay. 


r*  m 


In  addition  to  the  species  Eq.  (3.1)  or  its  numerical 
form  (3.5),  the  streamtube  calculation  requires  the  following 
equations 


du  m  _  1  d£ 
dt  p  dx 


h  -  H  -  \\*/2 


h  (T)  «  £  h£(T) 

i 


(3.9) 


p  -  p/tRTEOI^) 


(3.10) 


where  x  is  the  physical  distance  along  the  streamtube  and  all 
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differential  equations  have  ween  cast  in  terms  of  time,  t,  as 


the  independent  variable.  The  inputs  to  the  program  are  the 

initial  mass  fractions  a.,  velocity  u,  position  x,  and  either  the 

1 


, 


temperature  T  or  static  enthalpy.  Depending  upr i  what  is  ini- 

‘4 


tially  given,  Eq.  (3.9),  which  defines  the  static  enthalpy  in 


terms  of  the  mass  fractions  and  their  enthalpies,  is  solved  for 


the  initial  temperatures  or  static  enthalpy.  The  constant 


stagnation  enthalpy  H  is  then  determined  from  Eq.  (3.8)  and  the 


density  p  from  Eq.  (3.10).  The  pressure  terms  in  Eqs.  (3.6) 


and  (3.10)  are  given  input  functions  which  for  the  purpose  of 


this  program  are  assumed  to  be  quadratic  polynomials  in  x  over 


various  segments  of  the  streamtube.  Equation  (3.1)  is  then 


linearized  and  solved  for  y^  over  a  given  time  step.  If  the 
density  is  assumed  constant  over  a  step,  E-js.  (3.6)  and  (3.7) 


represent  a  linear  system  for  u  and  x,  which  is  also  solved  by  the 


Pade'  integration  procedure.  Using  the  new  values  of  y^  u,  and 


x,  the  thermodynamic  state  can  be  determined. 


An  Important  problem  in  any  integration  procedure  is 


the  selection  or  control  of  stepsize.  While  it  is  economically 


desirable  to  use  as  large  an  interval  as  possible,  the  stepsize 


must  not  be  increased  to  a  point  where  truncation  errors  may 


mask  the  true  solution.  Further,  the  stepsize  should  bo 
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the  chemical  system.  The  first  obvious  test  is  to  reduce  th* 
stepsize  wherever  negative  mass  fractions  appear.  However,  th 
alone  is  not  enough  since  the  mass  fraction  may  be  very 
inaccurate  although  not  negative.  The  criterion  in  the  stream- 
tube  program  is  to  determine  stepsize  on  the  basis  of  tempera¬ 
ture  change.  Experience  shows  that  in  regions  of  relatively 
small  temperature  change  the  solution  was  fairly  linear,  while 
it  .  strongly  non-linear  in  regions  of  great  temperature  change. 
Hence  temperature  change  could  be  used  as  a  means  of  controlling 
stepsize.  If  the  temperature  change  becomes  too  large  the 
stepsize  is  cut  in  half  and  if  the  temperature  change  is  small 
the  stepsize  is  doubled.  This  method  was  used  very  successfully 
with  dissociating  air  systems. 

Another  important  consideration  in  the  streamtube 
program  is  the  effect  of  retaining  only  first  order  species  terms 
in  the  numerical  method.  If  Eq.  (3.1)  is  written  in  the  general 
form  y  ■  f(y,T, p)  and  expanded  through  first  order  terms,  the 
linearized  equation  becomes 

y  -  f(y0'To'po)  +  (af/&y)0  Ay  +  (af/av)  at  +  (df/op)  a p 
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where  y,  f,  df/dT  and  df/&p  are  vectors,  &f/dy  is  ti  metrisi.  and 
(  )Q  implies  the  function  is  evaluated  at  the  beginning  of  the 
step.  The  density  term  can  oe  eliminated  in  terms  of  T  and  y  by 

introducing  the  equation  of  state.  However,  experience  has 

\  *  ' 

indicated  that  the  density  term  in  (3.11)  was  of  minor  impor- 
tance  in  the  numerical  solution  for  a  variety  of  chemical 
systems.  Hence  the  entire  Ap-term  is  eliminated  from  (3.11). 

The  method  of  handling  the  temperature  term  is  not.  straight¬ 
forward  and  several  procedures  have  been  utilized.  For  some 
chemistry  systems,  t  e  change  in  temperature  over  an  integration 
step  is  su  'ficiently  small,  and  the  AT-term  can  be  eliminated 
without  any  significant  effect  upon  the  numerical  solution. 
Equation  (3.11)  then  simplifies  to 

y  “  f(V0-T0-P0)  +  (jf)  (3.12) 

But  this  simplification  ie  only  satisfactory  for  some  air  chem¬ 
istry  systems  when  very  small  stepsizes  are  used.  1  overcome 
this  difficulty,  a  program  was  developed  to  include  the  term 
cdf /&T)o  AT  in  (3.11).  Initially  AT  is  unknown  and  the  program 
has  been  designed  to  iterate  until  tha  solution  converges  on 
the  correct  value.  As  mentioned  under  Subsection  3.1,  the 
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generator  program  haa  the  option  of  producing  programs  with  or 


without  the  added  temperature  cerm.  For  those  chemistry  systems 


which  do  not  involve  large  changes  in  temperature  over  an  inte¬ 


gration  step,  the  user  can  select  the  simpler  option  and  avoid 


the  added  computation  associated  with  iteration.  For  air 


chemistry  systems  as  illustrated  by  Appendix  B,  the  iteution 


option  should  be  included. 
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4.  results  from  streamtube  calculations 

While  the  technique*  of  numerical  analy.i.  have  been  uaed 
extensively  to  examine  truncation  error  and  obtain  indication* 
of  the  behavior  of  the  propagated  error,  it  is  not  possible  to 
adcurately  predict  the  error  in  a  numerical  solution  for  a 
given  problem.  In  some  cases,  a  bound  on  the  solution  error 
can  be  derived,  but  unfortunately  these  bounds  are  much  too 
conservative  to  provide  useful  information.  Consequently,  the 
solution  error,  which  is  an  important  consideration  in  the 
application  of  the  Pade  integration  procedure,  has  been  exam¬ 
ined  via  numerical  experimentation. 

Several  A  fluid  flow  problems  involving  various  chem- 

al  kinetics  systems  have  been  investigated  using  the  Pade* 
procedure,  and  the  results  have  been  compared  with  the  solution 
obtained  by  othei  numerical  methods.  Two  of  these  problems  are 
presented  herein.  For  both  problems,  the  Runge-Kutta  integra¬ 
tion  procedure  has  been  used  for  comparison  purposes.  Since  the 
Runge-Kutta  procedure  is  limited  with  respect  to  its  stepsize, 
the  numerical  solution  is  initially  obtained  using  a  stepsize 
which  is  considerably  less  than  predicted  by  the  stability 
analysis,  and  then  the  problem  is  rerun  using  a  stepsize  very 
close  to  the  stability  limit. 


In  this  manner  a  very  accurate 


solution  is  obtained  for  the  purposes  of  the  comparison  with  the 
results  from  the  Pade*  procedure.  The  second  run  of  the  Runge- 
Kutta  procedure  is  used  to  meauure  the  minimum  amount  of  machine 
time  to  obtain  the  solution.  In  this  way,  the  saving  in  machine 
time  through  the  use  of  the  Pade*  procedure  is  demonstrated 
without  any  built-in  penalty  against  the  Runge-Kutta  method. 

The  first  problem  to  be  considered  is  a  streamtube  calcula¬ 
tion  for  seven  species  air  consisting  of  0, ,  Na  ,  0,  N,  NO,  N0+, 
and  e~ .  The  initial  conditions  and  the  polynomials  describing 
the  variation  of  the  pressure  are  given  in  Table  I.  The  pressure 
polynomials  which  apply  over  discrete  intervals  are  expressed  in 
terms  of  s,  the  nondimensional  distance  along  the  streamtube. 

In  Fig.  2,  the  variation  of  temperature  and  veTjcity  along 
the  streamtube  is  shown.  The  sjlid  curves  «*i.c  results  obtained 
using  the  Runge-Kutta  integration  procedure  whereas  the  symboled 
points  are  from  the  Pade'  numerical  solution. 

In  Figs.  3  and  4,  the  behavior  of  the  species  a  and  a  _ 

U*j  6 

are  shown.  Initially  the  species  ctr  decreases  rapidly  from  its 
value  at  x  =  0.  In  contrast,  the  mass  fraction  of  electrons 
initially  increases  to  a  value  of  approximately  4.25  x  10  *  and 
then  starts  decreasing.  The  increase  in  electrons  occurs  over 
the  first  few  steps  and  represents  an  adjustment  in  the  initial 


46 


condition  (ae„  ^  0)  which  is  inconsistent  with  the  specified  temp¬ 


erature. 

In  all  of  these  figures,  the  results  from  the  Pade'  inte¬ 
gration  procedure  are  in  good  ugreement  with  the  results  from 
the  Runge-Kutta  procedure.  Regarding  the  amount  of  computer 
time,  the  Pade'  procedure  required  approximately  8  st»-  on  the 
IBM  7094  to  integrate  over  the  0.5  ft  domain,  whereas  the 
Runge-Kutta  method  took  200  sec.  Hence  for  this  chemistry  sys¬ 
tem  the  running  time  is  reduced  by  a  factor  of  about  25. 

As  a  second  example,  consider  the  chemistry  ryttem  for  the 
species  0,,  H,,  0,  N,  NO,  NO+,  «*,  0+,  0+ .  N+,  end  Na+.  Again  a 
streamtube  problem  was  run  using  the  Runge-Kutta  and  the  Pade' 
integration  procedures.  The  initial  conditions  and  the  pressure 
polynomials  are  given  xn  Table  I.T. 

In  Fig.  5,  the  temperature  and  velocity  results  are  pre¬ 
sented,  and  the  results  from  the  Runge-Kutta  (solid  line)  and 
the  Pade'  (symboled  points)  procedures  are  in  good  agreement. 
Similarly  the  species  results  from  the  two  numerical  procedures 
are  in  satisfactory  agreement  (Figs.  6  and  7). 

For  tne  mass  fraction  at  (Fig.  6)  the  scale  on  the  ordinate 
has  been  greatly  enlarged  to  illustrate  the  behavior  of  the 
Runge-Kutta  procedure  when  an  integration  step  is  sligntly  too 


I 


.  v , 


large.  The  value  of  the  mass  fraction  (at  x  -  .16)  has  an  oscix 


latory  kehavicr  and  If  the  stepslze  Is  not  Immediately  reduced 


the  amplitude  of  the  oscillation  would  rapidly  grow  in  an 


unbounded  manner.  At  such  a  point,  some  integration  routines 


cut  the  stepslze  and  restart  the  integration  procedure  at  an 


earlier  location  where  the  solution  was  stable.  For  the  purpose 


of  t'nis  study,  a  simpler  approach  was  used;  namely,  the  inte¬ 


gration  routine  detects  the  instability,  reduces  the  stepslze, 


and  continues  the  integration. 


For  this  problem,  it  was  convenient  to  make  the  comparison 


for  computer  time  over  the  streamtube  interval  (0,  0.2)  ft.  The 


Pade’  procedure  required  28  sec  whereas  the  Runge-Kutta  about 


650  sec.  Extrapolating  these  results  to  the  full  length  of  the 


streamtube  calculation  (.5  ft),  the  running  times  would  be 


70  sec  and  25  min  for  the  Pade'  and  Runge-Kutta  methods,  respec¬ 


tively.  The  example  serves  to  illustrate  the  importance  of  the 


rational  approximat.  on  method  in  reducing  the  amount  of  computer 


time  for  streamtube  calculations. 
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TABLE  I 


INITIAL  CONDITIONS  AND  PRESSURE  VARIATION 
FOR  THE  SEVEN  SPECJES  AIR  PROBLEM 


T  *  16836  0  K 

U  -  6524  ft/sec 


r'-.rjti 


.2378 

.7622 


“n  "  «NO  "  W  “  V  "  0 


Pressure  Polynomials 


p  -  aQ  +  a;  s  +  a,s8  (lb/ft*) 


where  s  ■  x/t,  x  is  the  distance  along  the  streamline  and 
l  *  0.5  ft,  the  nose  radius. 


ao 

as 

Domain  of  s  for 

Start 

Polynomial  Fit 

End 

5961.1 

-3694.7 

0 

0 

0.19059 

6817.1 

-8186.1 

0 

0.19059 

0.551 

6181.1 

-9100.2 

3753.9 

0.551 

0 . 681 

6344.7 

-9318.3 

3/21.4 

0.681 

0.8168 

6102.0 

-8699.7 

3326.6 

0.8168 

1.0 

..  > 
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TABLE  II 


INITIAL  CONDITIONS  AND  PRESSURE  VARIATION 
FOR  THE  ELEVEN  SPECIES  AIR  PROBLEM 


T  -  22806  K 

U  -  7561  ft/aec 


%  "  *232 


-  .768 


x  +  ■  a  ■  a  , 
NO+  •-  o+ 


Pressure  Poj.-r..jmiale 
(See  Table  I  for  definition  of  Polynomial  form) 


Domain  of  ■  for  Polynomial  Fit 


Start 


9671.3 

-  5994.3 

0 

0 

0.19059 

• 

11060.1 

-13280.2 

0 

0.1905,) 

0.551 

10028.2 

-14764.1 

8090.4 

0.551 

0.681 

V 

■  , 

10293.7 

-15118.1 

6037.6 

0.681 

0.8168 

►  **  ’«  j  •  ;  * 

*  " 
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9901.2 

-14114.4 

5397.1 

0.8168 

1.0 
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It  was  shown  in  Section  2  that  standard  methods  for  inte¬ 


grating  ODE  place  a  restriction  on  the  stepsize.  Unless  the  step- 
size  is  kept  relatively  small,  the  numerical  solution  will  become 
unstable.  As  shown  below,  similar  instabilities  arise  in  the 
numerical  solution  of  partial  differential  equations.  However, 
a  much  less  severe  restriction  upon  stepsize  can  be  obtained  by 
suitably  transforming  the  equations  and  then  applying  the  methods 
of  Section  2. 

For  standard  finite  difference  methods  the  numerical  solu¬ 
tion  of  a  parabolic  partial  differential  equation  will  not  be 
stable  unless  certain  restrictions  are  placed  on  the  mesh  size. 
Consider  an  equation  of  the  form: 

afc(x,t)  -  a  axx(x,t)  -  b  a(x,t) 
a(x,o)  -  f(x) 


(5.1) 


a,b  >  0 

where  f  (x)  and  the  solution  Ot(x,t)  are  bounded.  Let  the  numerical 
solution  be  denoted  by  A(x,t).  Place  a  uniform  .tiesh  size  of  h  in 
the  x  direction  and  K  in  the  t  direction  with  lattice  points 
(xi#  t^)  where  xi  -  ih;  i  -  0,  ±1,  ±2  ...  ,  and  t.,  »  jK;  j  -  0, 

1,  2  ...  Further  set  A(x  t.)  “A,  ,.  Using  forward  and  central 
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difference*  the  numerical  idution  mu*t  satisfy  the  finite  dlffe  - 
«n«  equation 


(Ai,j+i  -  Vj'  <*1+1,1  - 2  *1, 


(5.2) 


Ai,j+1  *  X*  Ai+l,j  *  (1  -  2  X»  -  2<b)A^  +  la  A 


i-l.J 


where  X  -  K/h* .  Starting  with  A£  -  (5.2)  can  be  solved 


for  A 


i,j+l  knowin9  Ai/ j*  However,  if  Kb  >  2  then  since  a  >  0 


(5.3) 


and  an  error  introduced  into  the  numerical  solution  will  not 


remain  bounded.  To  illustrate  the  numerical  behavior  let 

f 

«?  i  ■  0 

A.  ■  - 

0;  i  fi  0 


The  solution  of  (5.2)  corresponding  to  this  initial  condition  will 
show  how  a  single  error  introduced  at  i*0  propagates.  First  note 
that  if  A  is  non-zero  at  any  two  adjacent  points  it  differs  in 
sign  at  these  points.  This  is  obvious  for  all  i  at  j-0  and  j-i. 
From  (5.2)  and  (5.3)  it  follows  immediately  that  if  A.  has  one 

*■*  j 

sign  then  A1  H+1  will  have  the  opposite  sign.  Hence 


''  >  ■  • 

•S'  *  •  %  J 


/ 


j 
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For  any  j  at  moat  2j+l  terras  are  non-zero.  Suppose  the  term  of 
largest  magnitude  is  A^.  Then  *  *2j  +  1)  |An  |  and 

|An  j  |  >  (4Xa  +  103  “  */ (2 j  +  1).  But,  if  Kb  >  2, 

4Xa  +  Kb  -  1  >  1.  Hence,  the  initial  error,  is  amplified  from 
step  to  step  and  the  solution  is  unbounded.  For  the  numerical 
method  to  remain  stable,  the  necessary  condition  is 


or 


1  -  4Xa  -  Kb 


K  <  — Zh! 

K  4a+h*b 


<  1 


If  hab  «  a,  the  well  known  stability  condition  for  the  simple 
diffusion  problem  follows  from  (5.4).  However,  in  the  usual 
finite  rate  problem,  the  b  terra,  which  represents  the  species 
generation,  is  very  large  compared  tc  a.  Consequently,  the 


stability  condition  ia  approximately 


' .  7  ; 

•  f 


K  <f 


(5.5) 


This  limitation  upon  the  stepsize  ia  the  aame  form  as  discussed 


in  faction  2  for  ordinary  differential  equationa  uaing  the 


Euler  method.  For  a  large  value  of  h,  the  stepaize  ia  aeverely 


reatricted. 


The  above  result  has  been  derived  using  an  explicit  type  of 


difference  formulation  for  (5.1),  and  it  in  quite  natural  to 


try  and  overcome  the  restriction  with  an  implicit  differencing 


'  a*  ~ 

*  I 


method.  For  the  simple  diffusion  problem  (b^O),  the  method  ia 


known  to  be  unconditionally  stable  and  by  a  direct  extension. 


the  method  can  be  shown  to  be  stable  when  b  /  0.  However,  the 


simple  form  of  Eq.  (5.1)  is  not  representative  of  species  con¬ 


servation  equations  when  implicit  numerical  methods  are  used. 


The  species  generation  terms,  which  have  the  form  of  the  right 


side  of  Eq.  (3.1),  are  non-linear  functions  of  all  the  dependent 


variables.  Within  a  suitable  approximation,  these  terms  can  be 


r-, 


linearized  as  described  in  Section  2,  and  then  each  species 


conservation  equation  is  linear  but  still  coupled  with  every 


other  species.  Specifically,  (5.1)  becomes 


m 


#  • 

*  i 


i 

‘i. 
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a{tl)  (x.t) 


a  (x,t) 

xx 


n 


-  I 


b(k,«(k)(*.t) 


-  c 


(r) 


k-1 


(5.6) 


where  the  superscript  r  is  an  index  ranging  over  the  number  of 

(]^  \  ( r ) 

species  and  the  coefficients  b  and  c  are  obtained  from  the 
linearized  i'orm  of  (3,1).  The  implicit  method  of  solution  of 
(5.6)  can  be  applied  in  several  ways.  For  example,  if  the 
integration  is  from  t_.  to  t^+^  >  the  difference  scheme  can  be 


written  as 


aX  +  ( 1  +  b  ^  K  +  2al)A^r^  -  aX  A^r^ 

A  i-l,j+l  1  K  2  A)Ai,j+l  **  Ai+l,j+l 


*  5 

k?i 


Kb(k)A.<k>  , 

1#  j+1 


A<r>  -  Kc(r) 

1*  j 


(5.7) 


If  there  are  n  species  and  i  lattice  points  in  the  x  direction 
the  system  of  simultaneous  equations  involves  n*i  unknowns.  The 
coefficient  matrix  is  tridiagonal  along  the  main  diagonal  with 
n-1  additional  diagonals  spaced  every  i  columns.  Hence  the 
matrix  has  numerous  zero  elements.  Although  algorithms  could  be 
developed  to  efficiently  solve  the  systems  of  equatior ,  the  number 
of  machine  operations  would  still  be  large,  and  since  the  calcu¬ 
lations  must  be  repeated  at  each  integration  step,  the  amount 
of  machine  time  may  become  excessive.  Furthermore,  the 
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amplification  matrix,  which  is  the  inverse  of  the  matrix  repre¬ 
sented  by  the  left  side  of  (5.7),  miay  only  be  convergent  with 

a  small  stepsize.  Unlike  the  usual  diffusion  problem,  the  off- 

(k) 

diagonal  terms  b  (Mj)  can  be  very  large,  and  the  stability 

behavior  of  the  implicit  procedure  is  cor.tr letely  titered.  There 

fore,  this  numerical  procedure  is  not  considered  satisfactory 

for  the  efficient  integration  of  the  reaction  equations.  This 

conclusion  has  also  been  verified  by  numerical  experimentation 

with  small  chemical  kinetic  systems. 

A  second  method  of  implementing  the  implicit  method  is  to 

(k) 

iterate  for  the  terms  involving  the  b  coefficients.  The 
difference  equation  becomes: 


-  aX  A 


(r) 

i-1*  j+1 


(1  +  Kb^ 


+ 


2Xa)A 


(r) 

i#  j+1 


-  aX  A 


(r) 

i+1#  j+1 


-  A(r) 

i,j 


-  Kc 


(r) 


-  z 


Mr 


(k)  (k) 

Kb 


(5.8) 


and  the  matrix  is  then  tridiagonal  and  easily  inverted.  The 
A^  terms  on  the  right  side  are  initially  evaluated  using  the 
numerical  results  from  the  previous  station;  i.e.  A^.  Subsequent 
iterations  utilize  the  previouw~y  iterated  value  of  A.  .  . . 

1  #  J 

Unfortunately,  this  iteration  procedure  is  convergent  only,  if 
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the  stepsize  K  is  small;  again,  the  cause  of  this  difficulty  is 
the  large  value  of  b^. 

To  overcome  the  numerical  difficulties  that  accompany  the 
standard  explicit  and  implicit  numerical  methods,  a  technique 
(dee  Ref.  6)  has  evolved  wherein  the  diffusion  and  species  gener¬ 
ation  terms  are  handled  independently.  The  generation  terms 
appear  as  ordinary  differential  equations,  and  the  diffusion 
» erms  *re  separate  partial  differential  equations.  The  solutions 
are  then  combined  at  the  end  of  each  integration  step.  Although 
the  procedure  imposes  rather  strong  assumptions  regarding  the 
species  behavior  over  a  step,  a  substantial  increase  in  the  step- 
size  is  possible.  For  this  method,  let  ot  *  p  +  q,  where  p  and 
q  are  defined  between  t^  and  t^+1  by 

Pt(x,tj  -  a  «*xx(x,t)  (5.8) 

qt(x,t)  -  -  b  <*(x,t)  (5.9) 

In  order  to  solve  (5.9)  at  x  ■  assume  the  p  term  is  constant 
between  ^  and  t j+1;  that  is,  pfx^t)  *  p.  .  Then  in  (5.9) 
the  only  independent  variable  is  t,  and 


t 


Jl  ,  A  i 
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q'(t)  -  -  b  qi(t)  -  b  Pj( 

qi(°)  -qM  -0 


(5.10) 


The  solution  of  (5.10)  is 


-bt 


qiU)  "  ’  1} 


Equation  (5.8)  can  be  solved  directly  by  the  explicit  finite 
difference  method 


,j+:  "  Xa  Pi+i, 


j +  u  -  2Xa)pi.j +  Xa 


By  combining  these  two  solutions,  the  consistent  finite  difference 
form  is 

Ai,j+1  "  Xa  Ai+l,j  +  (®  “  2Xa)Ai,j  +  Xa  Ai_l,j 

Following  the  same  stability  argument  as  before,  if 


2Xa  >  1  >  e_Kb  -  e_Kb  -  2^  <  0 


and  for  every  j  there  is  a  point  (x  ,t.)  where 

n  j 


-Kb„ 


A  .  i  (4X«  -  e"103)?  «/ (2 j+1)  , 

n,  3 


Since  4Xa  -  e  >1,  an  error  committed  at  any  step  v/ill  grow. 
For  the  system  to  remain  stable 


'Vv>.. 
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or 


X  <  1/2 a 
K  <  h*/2a 


(5.12) 


This  criterion,  which  corresponds  to  the  standard  diffusion  con¬ 
dition,  is  much  less  restrictive  than  (5.5)  since  b  »  a.  Jhuu 
if  the  ordinary  differential  eq.  (5.9)  can  be  solved  accurately 
using  a  large  step,  the  numerical  method  can  overcome  the 
restriction  imposed  by  the  large  value  of  b.  While  standard 
methods  for  numerically  solving  the  ordinary  differential  eq.  (5.9) 
become  unstable  with  a  large  stepsize,  the  Pade'  integration 
method  developed  in  Section  2  will  yield  the  desired  results- 
The  prev/ous  stability  analysis  is  valid  with  e”103  replaced  by 
Q  1  (-Kb)  P (-Kb) .  Since 


<fl  (-Kb)  P  (-Kb) 


<  1 


the  overall  method  remains  stable  within  the  restriction  of 
(5.12) . 

The  integration  procedure  implied  by  Eq.  (5.11)  permits  a 
significantly  larger  stepsize  than  the  standard  explicit  or 
implicit  methods,  and  it  has  been  successfully  applied  to  several 
two-dimensional  problems.  The  saving  in  machine  time  for  these 


problem*  wa*  very  significant  and  further  saving*  are  expected 
a*  work  continues  in  the  study  and  application  of  th«.  method,  h 
problem  ol  particular  interest  is  the  relationship  between  the 
stable  stepsize  prescribed  by  (5.12)  and  the  limitation  imposed 
by  the  truncation  error.  Under  some  circumstances,  a  fraction 
of  the  stable  stepsize  must  be  o.i  d  to  achieve  the  desired 
accuracy  in  the  numerical  solution. 


™  L 


60 

6.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  limitations  on  the  stepsize  for  the  numerical  integration 
of  the  finite  rate  equations  has  been  investigated  and  expl  .ined 
in  terms  of  the  eigenvalues  of  the  ~mplif ication  matrix.  For 
standard  numerical  methods  of  integration  (Runge-Kuttc,  Predictor- 
Corrector,  etc.)  the  restriction  is  severe  and  causes  an  excas- 
sive  use  of  machine  time.  To  overcome  the  difficulty,  an  inte¬ 
gration  procedure  which  utilizes  a  rational  approximation  to 
replace  the  exponential  term  is  recommended.  This  method  is 
unconditionally  stable,  and  the  only  restriction  upon  the  step- 
size  is  imposed  by  the  usual  need  for  control  of  truncation  and 
round-off  errors  which  are  common  to  all  numerical  approximations. 

The  method  is  directly  applicable  to  linear  ordinary  differ¬ 
ential  equations  and  is  adapted  to  the  nonlinear  finite  rate 
equations  by  a  Taylor's  Series  expansion  in  n-variablea.  Since 
such  an  expansion  is  algebraically  difficult,  special  algorithms 
have  been  developed  and  programmed  to  automatically  perform  the 
expansion  and  produce  the  required  digital  computer  subroutines. 
The  generator  program  is  completely  operational  and  can  be  used  to 
develop  subroutines  for  generic  chemical  systems. 

The  numerical  procedure  has  been  used  in  the  analysis  of 
numerous  chemical  Kinetics  systems  including  systems  with  as  many 
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as  31  species  and  78  reactions.  For  an  air  chemistry  system  con¬ 
sisting  of  11  species,  the  rational  approx .imat ion  method  of  inte- 
9ra^on  approximately  1/20  the  machine  time  required  by  the 

Runge-Kutta  procedure.  Hence,  for  a  typical  stream  tube  problem, 
the  machine  time  was  reduced  from  25  minutes  to  70  seconds. 

The  method  has  been  investigated  for  application  to  partial 
ifferential  equations  which  can  not  be  satisfactorily  integ  ted 
using  either  the  standard  explicit  or  implicit  methods  because 
of  severe  stepsize  restrictions.  Two-dimensional  flow  problems 
involving  chemical  kinetics  are  of  this  type.  The  method  has  been 
applied  to  these  problems  and  satisfactory  results  have  been 
obtained.  However,  future  research  effort  should  be  directed 
toward  further  study  of  the  method  in  this  application.  In 
particular,  the  effect  and  control  of  truncation  error  should  be 
studied  in  detail  with  a  view  toward  developing  criteria  for  the 
control  of  the  integration  step. 


II 
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